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This  article  presents  hybrid,  degradation-based  reliability  models  for  a  single-unit  system  whose  degradation  is  driven  by  a  semi- 
Markov  environment.  The  primary  objective  is  to  develop  a  mathematical  framework  and  associated  computational  techniques  that 
unite  environmental  data  and  stochastic  failure  models  to  assess  the  current  or  future  health  of  the  system.  By  employing  phase-type 
distributions,  it  is  possible  to  construct  a  surrogate  environment  process  that  is  amenable  to  analysis  by  exact  Markovian  techniques 
to  obtain  reliability  estimates.  The  viability  of  the  proposed  approach  and  the  quality  of  the  approximations  are  demonstrated  in  two 
numerical  experiments.  The  numerical  results  indicate  that  remarkably  accurate  lifetime  distribution  and  moment  approximations 
are  attainable. 
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1.  Introduction 

Recent  advances  in  sensor  technologies  have  dramatically 
accelerated  the  use  of  sensors  for  continuously  monitor¬ 
ing  critical  engineering  components  in  a  host  of  applica¬ 
tions  (e.g.,  manufacturing  equipment,  aircraft  components, 
structures,  roadway  pavement,  power  grids,  etc.).  The  ad¬ 
vantage  of  using  sensor  data  to  estimate  the  current  or 
future  health  of  a  component  or  system  is  obvious — such 
data  obviate  the  need  for  failure  time  observations  which 
may  be  scarce.  For  example,  failure  time  data  may  be  very 
limited  when  systems  are  highly  reliable  or  when  it  is  pro¬ 
hibitively  expensive  to  run  systems  to  failure.  Consequently, 
many  researchers  have  been  advocating  degradation-based 
techniques  as  an  alternative  to  the  failure-based  paradigm 
in  order  to  exploit  a  plethora  of  degradation-related  data 
that  are  now  attainable  through  advanced  sensing  technolo¬ 
gies.  However,  in  some  applications,  the  degradation  may 
be  difficult  to  measure,  but  the  cumulative  degradation  can 
often  be  characterized  as  a  function  of  the  environment 
in  which  the  component  resides  or  the  conditions  under 
which  it  operates.  One  example  is  the  degradation  of  pro¬ 
tective  chemical  coatings  that  are  exposed  to  time-varying 
environmental  conditions.  If  the  environment  can  be  char¬ 
acterized  appropriately  (by  some  stochastic  process),  then 
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it  may  be  possible  to  characterize  the  cumulative  degrada¬ 
tion  as  a  function  of  time.  However,  translating  environ¬ 
mental  data  into  useful  lifetime  distributions  for  reliability 
estimation  remains  a  significant  challenge.  Therefore,  there 
exists  a  critical  need  for  novel  techniques  that  can  map 
environmental  sensor  data  to  degradation-based  reliability 
indices. 

In  this  article,  we  present  hybrid,  degradation-based  reli¬ 
ability  models  for  a  single-unit  system  whose  degradation  is 
driven  by  a  semi-Markov  environment.  The  primary  objec¬ 
tive  is  to  develop  a  mathematical  framework  and  associated 
computational  techniques  for  uniting  environmental  data 
and  stochastic  failure  models  to  assess  the  current  or  future 
health  of  the  system.  Specifically,  we  develop  a  general  pro¬ 
cedure  that  uses  information  about  the  sensed  environment 
to  estimate  important  reliability  indices  (e.g.,  the  reliabil¬ 
ity  function,  the  residual  lifetime  distribution,  the  mean 
time-to-failure  and  the  residual  mean  time-to-failure).  Our 
procedures  are  most  suitable  when  either;  (i)  discrete  (and 
distinct)  environment  states  can  be  identified  in  a  natural 
way;  or  (ii)  discrete  states  can  be  constructed  within  a  con¬ 
tinuous  state  space  in  a  natural  way.  The  framework  extends 
a  host  of  so-called  random  environment  models  to  the  case 
of  semi-Markovian  environments  that  place  only  mild  re¬ 
strictions  on  the  dynamics  of  the  evolving  environment. 
Moreover,  unlike  most  existing  stochastic  failure  models  of 
this  type,  we  devise  a  procedure  for  using  observed  data  to 
estimate  the  models’  parameters  and  to  produce  reliability 
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estimates.  The  need  to  estimate  reliability  indices  within  a 
non-failure-based  paradigm  has  become  critical  as  complex 
systems  are  increasingly  reliable  and  prohibitively  costly  to 
run  to  failure. 

The  literature  related  to  environment-  and  degradation- 
based  reliability,  while  relatively  immature  as  compared  to 
that  of  failure-based  reliability,  has  been  growing  at  a  rapid 
pace.  Generally,  work  in  this  area  can  be  partitioned  into 
two  categories  broadly  construed  as  physics-  and  statistics- 
based  models  and  probabilistic  (or  stochastic)  models. 
Physics-of-failure  models  are  typically  deterministic  mod¬ 
els  that  attempt  to  capture  the  basic  physical  principles 
underlying  the  failure  process.  Some  examples  are  pro¬ 
vided  in  Lawless  (1982),  Ebeling  (1997),  Mishra  and  Pecht 
(2002),  and  Elsayed  et  al.  (2006).  Statistics-based  mod¬ 
els  are  usually  developed  by  collecting  data  over  a  pe¬ 
riod  of  time  that  includes  a  reasonable  number  of  sys¬ 
tem  failures  and  then  using  these  data  to  form  the  basis 
of  a  statistical  model  that  predicts  failure  as  a  function 
of  the  sensor  measurement.  One  example  is  the  Propor¬ 
tional  Hazard  Model  (PHM)  developed  by  Cox  in  1972 
to  analyze  medical  survival  data  (Cox,  1972;  Cox  and 
Oakes,  1984).  The  PHM  has  been  implemented  in  a  variety 
of  engineering  applications,  such  as  aircraft,  marine  sys¬ 
tems,  and  machinery  (Jardine  and  Anderson,  1985;  Jardine 
et  al.,  1987,  1989;  Zhan  et  al.,  2003).  An  excellent  review 
of  PHMs  for  preventive  maintenance  can  be  found  in  Kob- 
bacy  et  al.  (1997).  Gebraeel,  et  al.  (2005)  and  Gebraeel 
and  Pan  (2008)  proposed  models  that  use  observed  sig¬ 
nals  to  fit  linear,  exponential,  and  polynomial  degradation 
curves  from  which  they  estimated  normal  residual  lifetime 
distributions.  The  results  were  compared  with  real  degra¬ 
dation  signals  using  frequency  measurements  from  rolling 
bearings.  Both  physics-of-failure  models  (by  definition)  and 
statistical  models  (due  to  the  experimental  methods  used  to 
develop  them)  are  limited  to  the  specihc  system  under  con¬ 
sideration.  Eurthermore,  because  it  may  be  difficult  and/or 
prohibitively  costly  to  run  systems  to  failure,  their  develop¬ 
ment  and  implementation  may  be  very  time-consuming,  as 
noted  by  Elsayed  (2000). 

Techniques  of  the  second  type  are  typically  based  on 
stochastic  shock  and  wear  models,  as  well  as  models  for 
systems  that  evolve  in  a  random  environment.  Early  work 
due  to  Esary  et  al.  (1973)  provided  several  results  for  both 
wear  and  shock  processes.  Qnlar  (1977,  1984)  generalized 
many  of  the  models  of  Esary  et  al.  (1973)  by  showing  that 
the  joint  process  of  the  unit’s  wear  level  and  the  state  of  its 
ambient  environment  may  be  viewed  as  a  Markov-additive 
process  and  provided  several  examples.  The  first  consid¬ 
ered  the  case  when  the  random  environment  is  a  finite 
Markov  process,  and  the  wear  is  assumed  to  increase  as 
a  Eevy  process.  Random  shocks  were  also  assumed  to  oc¬ 
cur  at  environment  transition  epochs.  The  second,  which  is 
similar  to  our  approach,  views  cumulative  wear  as  a  con¬ 
tinuous,  additive  functional  of  the  operating  environment, 
and  the  first  time-to-failure  is  a  first  passage  time  for  the 


cumulative  wear  process.  Other  stochastic  failure  models 
also  attempt  to  capture  the  impact  of  a  randomly  varying 
environment.  An  excellent  survey  of  such  models  is  due 
to  Singpurwalla  (1995).  Ei  and  Euo  (2005)  considered  a 
Markov-modulated  shock  process  wherein  the  shock  inter¬ 
arrival  times  and  the  random  shock  damage  are  both  gov¬ 
erned  by  a  Markov  chain.  They  obtained  reliability  bounds 
for  when  the  inter-arrival  times  have  heavy-  or  light -tailed 
distributions,  but  their  degradation  model  does  not  in¬ 
clude  a  continuous  wear  component.  Mallor  and  Omey 
(2001)  considered  a  generalized  shock  process  and  stud¬ 
ied  some  asymptotic  properties.  Kharoufeh  (2003)  and 
Kharoufeh  and  Cox  (2005)  presented  degradation  models 
that  assume  a  Markovian  environment.  Klutke  et  al.  (1996) 
examined  the  availability  of  an  inspected  system  whose 
inter-inspection  times  and  wear  rates  are  random.  Kiessler 
et  al.  (2002)  investigated  the  limiting  average  availability  of 
a  system  whose  time-varying  wear  rates  are  governed  by 
a  continuous-time  Markov  chain.  Kharoufeh  et  al.  (2006) 
extended  the  model  of  Kiessler  et  al.  (2002)  by  including 
damage-inducing  shocks  that  arrive  according  to  a  time- 
homogeneous  Poisson  process  and  deriving  the  Eaplace- 
Stieltjes  transforms  of  transient  reliability  indices.  Ebrahimi 
(2006)  considered  a  system  whose  degradation  is  comprised 
of  a  continuous  wear  component  as  well  as  jumps.  The 
properties  of  the  model  were  investigated  and  bounds  were 
established  for  the  reliability  function.  Recently,  Ozekici 
and  Soyer  (2003, 2004, 2006)  analytically  examined  reliabil¬ 
ity  indices  in  both  Markov  and  semi-Markov  environments. 

The  present  article  extends  the  models  of  Kiessler  et  al. 
(2002),  Gebraeel  et  al.  (2005),  Kharoufeh  and  Cox  (2005), 
and  Gebraeel  and  Pan  (2008)  in  very  important  ways. 
Eirst,  it  incorporates  environmental  data  for  the  purpose 
of  evaluating  reliability  indices  by  linking  these  data  to 
degradation.  Second,  it  generalizes  the  models  of  Kiessler 
et  al.  (2002),  Kharoufeh  (2003),  and  Kharoufeh  and  Cox 
(2005)  by  assuming  that  the  environment  does  not  evolve 
as  a  Markov  process  but  rather  as  a  semi-Markov  process. 
This  generalization  allows  us  to  consider  environment 
dynamics  that  are  not  restricted  to  memoryless  holding 
times  in  individual  environment  states.  Third,  it  provides 
a  viable,  data-driven  approach  for  translating  environment 
observations  into  failure  time  estimates  using  analytical 
stochastic  failure  models.  Specifically,  we  exploit  phase- 
type  (PH)  distributions  to  approximate  environment  state 
holding  times,  thereby  inducing  a  Markovian  structure 
that  is  amenable  to  exact  analysis  by  the  methods  described 
in  Kharoufeh  and  Cox  (2005)  and  Kharoufeh  et  al.  (2006). 
We  provide  guidance  on  selecting  an  appropriate  PH 
distribution  approximation  for  each  environment  holding 
time.  Einally,  using  simulated  benchmarks,  we  illustrate 
the  viability  of  our  approach,  and  the  remarkably  high 
quality  of  our  approximations,  by  way  of  two  numerical 
experiments.  The  key  innovation  of  our  work  is  the  use 
of  environmental  data  to  compute  lifetime  distributions 
within  a  degradation-based  paradigm. 
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The  remainder  of  the  article  is  organized  as  follows. 
Section  2  describes  the  semi-Markov  environment  model 
and  the  stochastic  failure  models.  In  Section  3,  we 
briefly  review  PH  distributions  and  provide  decision  rules 
for  selecting  a  PH  approximation.  Section  4  formally 
describes  our  procedure  for  converting  the  non-Markovian 
environment  into  one  that  is  amenable  to  exact  analysis  by 
Markovian  techniques.  Section  5  presents  two  randomized 
experiments  that  illustrate  the  viability  of  our  approach, 
while  Section  6  provides  some  concluding  remarks  and 
directions  for  future  work. 


2.  Model  description 

This  section  describes  our  mathematical  model  of  a  single¬ 
unit  system  (i.e.,  a  component)  subject  to  continuous, 
environment -driven  degradation.  The  component  is  placed 
into  service  at  time  0  in  perfect  working  order  and  degrades 
at  a  rate  that  depends  on  the  state  of  its  environment,  which 
evolves  as  a  continuous-time  stochastic  process  with  finitely 
many  discrete  states.  When  the  component’s  cumulative 
degradation  reaches  a  deterministic,  critical  threshold  value 
X,  the  system  is  said  to  be  failed.  We  denote  the  random  time 
to  first  reach  X  by  7^,  and  we  assume  that  is  proper  (i.e.,  as 
t  ^  oo,  P(Tx  <  t)  ^  1  for  each  x  >  0).  The  time-varying 
degradation  rate  is  modulated  by  an  external  stochastic 
process  commonly  termed  the  random  environment. 

The  random  environment  is  denoted  by  Z  =  { Z{t)  :  t  > 
0}  with  finite  state  space  5  =  { 1 ,  2, . . . ,  W}  and  2  <  K  <  oo. 
The  process  visits  some  state  i  e  S  and  spends  a  random 
amount  of  time  there  that  depends  on  the  next  state  it  will 
visit,  say,  j  e  S,  j  ^  i.  It  chooses  j  according  to  a  Markov 
chain  with  transition  probability  matrix  P.  The  time  spent 
in  state  i,  given  that  the  process  next  visits  state  j,  has  cu¬ 
mulative  distribution  function  (c.d.f )  Hij .  Let  denote  the 
time  of  the  «th  transition  of  the  environment  process,  and 
let  En  =  Z(5j+)  be  the  state  of  Zjust  after  the  «th  transi¬ 
tion  epoch.  The  environment  is  a  temporally  homogeneous 
Semi-Markov  process  (SMP)  with  associated  Markov  re¬ 
newal  sequence  [{E„  ,Sn)\n>  0}  on  the  state  space  S.  For 
each  t  >  0,  «  >  0  and  i,  j  e  S,  define  the  probabilities: 

Gijit)  =  P{En+\  =  j,  S„+1  -  Sn<t\En  =  i), 

and  the  corresponding  semi-Markov  kernel  matrix  G(?)  = 
[G, ,/(?)].  The  environment  transitions  from  state  i  to  state 
j  7^  i  according  to  a  Markov  chain  on  S  with  transition 
probability  matrix  P  where  the  (z,  yjth  element  of  P  is  given 
by 

Pij  =  P(Ei  =  j\Eo  =  i)  =  lira  (1) 

Although  the  transition  probabilities  {{/>,,/})  can  be  ob¬ 
tained  by  evaluating  the  limit  of  Equation  (1),  the  ker¬ 
nel  functions  Gij(l)  may  not  be  known  in  practice,  or 
they  may  be  difficult  to  characterize.  We  circumvent  this 


shortcoming  by  estimating  pi  j  from  observable  data  in 
Section  4. 

When  Z{t)  =  i  €  S,  the  c.d.f  of  the  environment’s  hold¬ 
ing  time  in  state  z,  given  that  it  next  visits  state  j  ^  i,  is 

=  P{Sn+i  -  Sn<t  \  En+\  =  j,  E„  =  z) 

=  P(5i  <  z  I  Fi  =  y,  Fo  =  i\ 

where  the  second  equality  implies  temporal  homogeneity. 
Therefore,  the  c.d.f  of  the  holding  time  in  state  z,  indepen¬ 
dent  of  the  next  state,  is 

mt)  =  P{S,  <  z  I  Fo  =  0  =  I]  Gijlt),  i  e  5. 

j€S 

For  Z  >  0  and  i,  j  e  S,  define  the  transition  functions  of 
Zby 

Jtijit)  ^  P(Z(t)  =  j  I  Z(0)  =  z), 

and  the  transition  matrix  P(z)  =  [7r,_y(z)].  It  can  be  shown 
(cf  Kulkarni  (1995))  that  7r,,y(z)  satisfies  the  Markov  re¬ 
newal  equation: 

7r,',y(z)  =  77, (z)  lj,=yj  +  V  /  nk,j(t  -  v)  Gi^ddv),  (2) 

where  77, (z)  =  1  -  77(z),  the  indicator  function  lj,=/)  as¬ 
sumes  a  value  of  one  if  z  =  j  and  is  zero  otherwise,  and 
the  integral  on  the  right-hand  side  of  Equation  (2)  is  the 
convolution  of  Jtkj  with  G/y,.  In  general,  it  is  difficult  to  ob¬ 
tain  the  transition  functions,  even  when  the  kernel  matrix 
is  known. 

Next,  define  a  rate  function  z-  :  5  ^  (0,  oo)  such  that 
whenever  Z{z)  =  z,  the  system  degrades  at  a  constant  rate 
r(i),  r{i)  >  0  for  z  =  1,  2, . . . ,  W.  That  is,  the  component’s 
degradation  rate  is  semi-Markov  modulated.  We  pause 
here  to  make  two  remarks.  First,  it  is  important  to  note  that 
our  technique  does  not  require,  or  assume,  that  the  overall 
degradation  path  is  linear  in  form.  Rather,  the  imposed 
assumption  is  that  the  rate  of  degradation  is  constant 
within  a  given  environment  state.  This  framework  allows  us 
to  characterize  the  degradation  process  by  a  set  of  mean 
rates  of  degradation  as  a  function  of  time  (or  usage). 
The  environment-modulated  rates  lead  to  piece-wise 
linear  sample  paths  of  degradation  that  can  approximate 
common  degradation  patterns  (e.g.,  linear,  exponential, 
polynomial  or  other  forms).  For  example,  once  a  crack 
has  been  initiated  in  a  metallic  specimen,  the  crack  grows 
exponentially  in  the  number  of  load  cycles  (see  Virkler  et 
al.  (1979)).  The  degradation  pattern  can  generally  be  rep¬ 
resented  by  three  distinct  regions.  The  first  region  exhibits 
mild  linear  growth  pattern  until  reaching  an  inflection 
point  at  which  the  degradation  grows  at  a  moderately 
higher  rate.  After  growing  moderately  for  some  time,  the 
degradation  path  reaches  a  second  inflection  point  beyond 
which  the  crack  growth  rate  is  very  high  until  reaching 
a  critical  failure  threshold  (e.g.,  the  fracture  point  of  the 
material).  It  was  shown  in  Kharoufeh  and  Cox  (2005)  that 
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such  a  nonlinear  degradation  path  can  be  approximated 
using  a  piece -wise  linear  path  with  three  distinct  rates 
associated  to  the  three  phases  of  crack  length  growth.  They 
developed  a  statistical  procedure  to:  (i)  estimate  the  appro¬ 
priate  number  of  environment  states;  and  (ii)  estimate  the 
degradation  rate  in  each  environment  state.  In  practice,  the 
degradation  rate  associated  with  each  environment  state 
must  be  estimated  with  great  care — in  consultation  with 
subject  matter  experts — to  ensure  consistency  with  phys¬ 
ical  principles.  Second,  the  assumption  of  strictly  positive 
degradation  rates  can  be  justified  as  follows.  It  is  presumed 
that,  whenever  the  unit  is  in  operation,  degradation  accrues 
to  at  least  some  degree;  i.e.,  the  component  experiences 
day-to-day  wear  from  normal  usage.  On  the  other  hand, 
in  some  scenarios,  even  if  the  system  is  not  in  use,  the 
ambient  environment  might  induce  degradation  (e.g., 
the  degradation  of  protective  chemical  coatings  that  are 
exposed  to  the  elements).  We  do  not  consider  environment 
states  that  improve  the  condition  of  the  component. 

The  cumulative  degradation  up  to  time  t,  denoted  by 
X{t),  is 

X(t)  =  X(0)+f  r[Z(u)]du,  (3) 

Jo 

where  we  assume  X(0)  s  0,  and 

/  \r[Z{u)]\du  <  oo  a.s., 

Jo 

so  that  X(l)  is  well  dehned  for  each  ?  >  0.  The  pro¬ 
cess,  X  =  {X{t) :  ?  >  0},  is  termed  the  cumulative  degra¬ 
dation  process.  The  positivity  of  the  degradation  rates, 
r(l),  r{2),  . . .  ,r(K),  ensures  that  the  sample  paths  of  X are 
almost  surely  monotone  increasing,  and  consequently,  that 
events  [X{t)  <  x}  and  {Ti  >  ?}  are  equivalent.  The  system’s 
random  lifetime  is  given  by 

=  inf{?  >  0  :  X{t)  >  x},  (4) 

or  the  first  time  the  degradation  process  X  reaches  x.  Let 

Fix,  t)  =  PiT^  <0  =  1-  PiMt)  <  X), 

denote  the  unconditional  c.d.f  of  the  unit’s  lifetime,  and  let 

7j(x,  t)=P{T^  <  1 1  Z(0)  =  0=1-  P{Xit)  <  X I  Z{0)  =  i), 

be  the  conditional  c.d.f  of  7^,  given  the  initial  state  of  the 
environment.  Denote  by  E[T^]  the  nth  moment  of  7^,  for 
n  >  1,  and  let  its  conditional  counterpart  be  denoted  by 
s  i?[r; I z(0)  =  /]. 

For  a  semi-Markov  environment  process  Z,  comput¬ 
ing  F{x,  t)  and  7j(x,  t)  is  non-trivial  due  to  the  difficulty 
in  obtaining  7r,_y(0  of  Equation  (2).  However,  if  Z  is  a 
Continuous-Time  Markov  Chain  (CTMC)  on  the  state 
space  S,  the  Laplace-Stieltjes  Transforms  (LSTs)  of  Fix,  t), 
Fiix,  t),  E[T^]  and  Ei[T^],  with  respect  to  x,  can  be  derived 
explicitly  (see  Kharoufeh  and  Cox  (2005)  and  Kharoufeh 
et  al.  (2006)).  Let  a  be  the  initial  distribution  vector  of  the 


environment,  e  is  a  column  vector  of  ones,  and  e,-  is  a  col¬ 
umn  vector  whose  ith  entry  is  unity  and  all  others  are  zero. 
The  unconditional,  full  lifetime  distribution  is  given  by 

Fix,  t)  =  PiTx  <0  =  1-  «V(x,  t)e 

where  V(x,  t)  =  [ ix,  t)]  is  a  K  x  K  matrix  with 

Vi,j  ix,  t)  =  PiXit)  <  X,  Zit)  =  7jZ(0)  =  /), 

the  joint  probability  that,  at  time  t,  the  degradation  of  the 
system  has  not  exceeded  x,  and  the  environment  process  is 
in  state  j  e  S,  given  that  the  environment  was  initially  in 
state  i  e  S.  The  conditional  c.d.f  of  7i  is 

Fiix,  t)  =  PiiTx  <0  =  1-  e;V(x,  Oe,  (5) 

where  e'  denotes  the  transpose  of  e, .  Let  the  LSTs  of  Fix,  t) 
and  Fiix,  t),  with  respect  to  x,  be 

POO 

Fiu,  t)  =  I  e~“^Fidx,  t),  Re(M)  >  0, 

Jo 

and 

POO 

Fiiu,  t)  =  I  e~“^Fiidx,  t),  Re(M)  >  0, 

Jo 

respectively.  Using  Theorem  3  of  Kharoufeh  et  al.  (2006), 
it  can  be  shown  that 

Fiu,  0  =  1-  aexp[(Q  -  wROOe,  (6) 

and 

Fiiu,  0  =  1-  ej  exp[(Q  -  MRd)f]e  (7) 

where  Q  is  the  generator  matrix  of  the  CTMC,  Rd  = 
diag(r(l),  r(2), . . . ,  r(K)),  and  exp(A)  denotes  exponenti¬ 
ation  of  the  square  matrix  A. 

While  these  results  are  useful,  they  are  applicable  only 
when  the  environment  evolves  as  a  CTMC  (i.e.,  when  the 
environment  spends  an  exponential  time  in  each  state  be¬ 
fore  it  transitions  to  another  state).  However,  in  reality, 
these  state  holding  times  are  general,  non-negative  random 
variables  whose  distributions  are  (in  most  cases)  unknown. 
In  fact,  the  assumption  of  memoryless  holding  time  dis¬ 
tributions  may  be  completely  unwarranted  in  many  appli¬ 
cations.  Our  aim  is  to  provide  a  means  by  which  to  use 
sensor  data,  especially  environment-related  data,  to  char¬ 
acterize  the  holding  time  distributions  and  the  environ¬ 
ment’s  transition  rates  so  that  Equation  (6)  or  Equation 
(7)  can  be  directly  applied  to  compute  reliability  indices. 
Our  technique  requires  observations  of  the  environment’s 
current  state,  its  transition  times,  and  observation  of  its 
subsequent  state.  These  will  be  used  in  an  automated 
procedure  to  estimate  holding  time  distributions  by  PH 
distributions.  By  doing  so,  we  generalize  the  models  pre¬ 
sented  in  Kiessler  et  al.  (2002),  Kharoufeh  (2003),  Gebraeel 
et  al.  (2005),  Kharoufeh  and  Cox  (2005),  and  Gebraeel  and 
Pan  (2008)  to  account  for  semi-Markovian  environments. 
Sections  3  and  4  describe  our  approach  for  approximating 
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the  c.d.f.  F{x,  t)  (or  ?))  and  its  moments  within  the  environment  that  is  Markovian  so  that  Equations  (6)  and 
semi-Markov  framework  presented  in  this  section.  (7)  can  be  used  to  compute  reliability  indices. 


3.  PH  approximations 

Here,  we  review  the  most  commonly  used  techniques  for 
approximating  the  c.d.f  of  a  non-negative  random  vari¬ 
able  by  a  PH  distribution.  For  a  complete  introduction  to 
PH  distributions,  the  reader  is  referred  to  Neuts  (1981). 
A  non-negative  random  variable,  Y,  is  said  to  have  a  PH 
distribution  if  it  represents  the  time  to  absorption  of  some 
Markov  chain.  If  the  chain  evolves  in  discrete  (continuous) 
time,  then  Eis  a  discrete  (continuous)  PH  random  variable. 
Owing  to  the  fact  that  our  environment  evolves  in  continu¬ 
ous  time,  we  review  only  the  continuous  version  here. 

Let  {(pit) :  ?  >  0}  be  a  CTMC  with  finite  state  space  M  = 
{1, 2, . . . ,  k  +  1}  where  states  1,  2, . . . ,  k  are  transient,  and 
state  k  +  1  is  absorbing.  The  infinitesimal  generator  matrix 
of  {(pit)  :  ?  >  0}  is 


where  T  is  a  A;  x  A;  matrix  with  2),,  <0  and  ./  >0,j^i, 
T°  is  a  column  vector,  0  is  the  zero  vector,  and  Te  +  T°  =  0. 
The  vector  (3o  =  (|3,  Pk+i)  is  the  initial  distribution  vector 
of  Q*;  i.e.,  |3  is  a  1  x  krow  vector,  and  Pt+i  is  the  probability 
that  the  Markov  process  begins  in  the  absorbing  state  k+  1. 
We  assume  throughout  that  |3o  =  0-0)-  The  vector  (3o 
satisfies  the  usual  condition  Pq  e  =  1 .  Let  Y  be  the  time 
to  absorption  of  Q* .  It  is  not  difficult  to  verify  (cf  Neuts 
(1981))  that  the  c.d.f  of  Tis 

Fiy)  =  P(  T  <  y)  =  1  -  p  exp(Ty)  e,  (9) 

and  the  nth  moment  of  Tis 


L[r]  =  (-l)"n!T^”e,  n  >  1. 

The  distribution  function  (9)  is  a  PH  distribution  with 
representation  (P,  T);  therefore,  constructing  a  PH  approx¬ 
imation  of  the  c.d.f  of  a  non-negative  random  variable 
X  involves  determining  the  pair  (p,  T).  PH  distributions 
are  attractive  because  they  can  approximate  the  c.d.f 
of  any  non-negative  random  variable.  Moreover,  they 
can  facilitate  tractability  when  analyzing  non-Markovian 
stochastic  processes.  That  is,  by  supplementing  the  original 
non-Markov  process  with  the  phase  variable,  (pit),  one  can 
obtain  a  Markovian  structure  so  that  standard  analysis 
techniques  for  Markov  processes  can  be  employed  (cf 
Whitt  (1984),  Altiok  (1985),  Perros  (1994),  and  Osogami 
and  Harchol-Balter  (2006)).  For  these  two  primary 
reasons,  we  will  employ  PH  distributions  to  approxi¬ 
mate  the  environment  state  holding  time  distributions, 
Fliit),  i  =  1,2, . . . ,  K.  By  doing  so,  we  will  replace  our 
semi-Markov  environment  by  a  surrogate  (expanded) 


3.1.  Specific  PH  Distributions 


The  general  k-phase  PH  distribution  (denoted  PHyt)  has 
an  associated  Markov  process  with  k  states  (or  phases) 
and  a  single  absorbing  state  k+l.  The  time  spent  in  a 
transient  state  i  is  exponentially  distributed  with  parameter 
IXi,i  =  1,  2, . . . ,  k.  Once  state  A:  +  1  is  reached  for  the  first 
time,  the  process  is  absorbed.  There  are  no  restrictions  on 
the  transitions  of  the  process  until  the  time  of  absorption. 
That  is,  a  visit  to  state  i,  I  <  i  <k,  can  be  followed  by  a 
visit  to  any  other  transient  state  before  absorption 

By  contrast,  the  A:-phase  Coxian  distribution  (denoted 
Ck)  is  the  distribution  of  the  time  to  absorption  of  a  Markov 
chain  with  k  phases  that  are  visited  sequentially,  but  ab¬ 
sorption  is  possible  at  the  end  of  each  phase.  As  for  the 
PHyt  distribution,  the  time  spent  in  phase  i  is  exponential 
with  rate  /x, ,  /  =  1 ,  2, . . . ,  Ar,  so  this  PH  distribution  has  the 
representation 


/— /xi  ai/xi  0  0 

I  0  — /X2  ^2/^2  0 


“  ) 
0 


T  = 


0  0 


0 


(10) 


Vo  0 


— M/t-1  ak-\tik~\  I 
0  -I^k  ) 


with  (3  =  (1,  0,  . . . ,  0)  and  T**  =  ((1  -  aO/xi,  (1  -  02)1^2, 
. . . ,  (1  —  ak-\)iak~\,  M/tX,  where  a,-  is  the  transition  prob¬ 
ability  from  A  to  A  +  1,  A  =  1,  2, . . . ,  A;  -  1,  and  1  -  a,  is  the 
probability  of  absorption  immediately  following  phase  A. 
The  Coxian  distribution  is  especially  attractive  because  it 
can  exactly  represent  any  distribution  that  has  a  rational 
LST  (Cox,  1955). 

The  Aj-phase  Erlang  distribution  (denoted  Ek)  is  a  spe¬ 
cial  case  of  the  Ck  distribution  in  which  /x,-  =  /x  for  A  = 
\,2, ..  .,k,  and  all  the  transient  phases  must  be  visited 
sequentially  before  being  absorbed  into  phase  k+  \  (i.e., 
Ui  =  1,  A  =  1,  2, . . . ,  A;  -  1)  so  that  absorption  is  not  possi¬ 
ble  until  after  the  Ath  phase  has  been  completed.  It  is  clear 
that  the  k-phase  Erlang  is  the  sum  of  k  independent  and 
identically  distributed  exponential  random  variables,  each 
having  parameter  /x.  Its  representation  (|3,  T)  is  given  by 


/— /X  /X  0  0  0  \ 

0  —/X  /X  0  0 

0  0  —/X  /X  0 


(11) 


Vo  0  0  ...  -/x/ 

with  p  =  (I,  0, . . . ,  0)  and  T**  =  (0,  0, . . . ,  0,  /x)'. 

The  Ar-phase  generalized  Erlang  distribution,  as  defined 
by  Marie  (1980),  is  identical  to  the  Ar-phase  Erlang  except 
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that  absorption  is  possible  at  the  end  of  phase  1.  That 
is,  {0{t)  :  ?  >  0}  transitions  from  phase  1  to  phase  2  with 
probability  a  and  is  absorbed  after  phase  1  with  probability 
\  -  a.  However,  if  phase  2  is  reached,  all  subsequent  phases 
are  visited  sequentially  until  absorption  into  state  k+\. 
Therefore,  the  representation  of  this  PH  distribution  is 


/— /X  a/x  0  0  0  ^ 

0  —/X  /X  0  0 

0  0  —/X  /X  0 


(12) 


Vo  0  0  ...  -/X/ 

with  (3  =  (1, 0, . . . ,  0)  and  T°  =  ((1  -  a)/x,  0, . . . ,  0,  /x)'. 

Any  one  of  the  PH/t,  Ck,  Ek,  and  ^-phase  generalized 
Erlang  distributions  can  be  used  to  approximate  the  c.d.f 
of  an  arbitrary,  non-negative  random  variable.  However, 
selecting  the  best  distribution  for  a  given  application  is 
not  straightforward.  Next,  we  provide  some  guidelines  for 
choosing  an  appropriate  approximation. 


3.2.  PH  approximation  selection 

PH  distribution  selection  is  a  critical  aspect  of  our  proposed 
technique,  and  three  important  criteria  must  be  considered 
when  making  the  selection.  First,  the  approximations  must 
use  the  least  possible  number  of  phases  while  adequately 
representing  the  environment  sojourn  times.  Second,  the 
approximation  should  accommodate  an  automated  proce¬ 
dure  so  that  only  observed  sensor  data  is  needed  to  con¬ 
struct  the  approximation.  Finally,  the  computational  effort 
should  be  minimal. 

Suppose  that  A  is  a  non-negative  random  variable  with 
unknown  c.d.f  H.  We  want  to  approximate  A  by  a  PH 
random  variable,  say  Y,  with  c.d.f  H.  Marie  (1980),  Altiok 
(1985),  Johnson  (1993),  Perros  (1994),  and  Osogami  and 
Harchol-Balter  (2006)  all  provide  excellent  summaries  of 
techniques  that  can  be  used  to  obtain  PH  distribution  ap¬ 
proximations.  With  the  exception  of  Osogami  and  Harchol- 
Balter  (2006),  most  techniques  require  an  estimate  of  the 
squared  coefficient  of  variation  of  X,  lm\,  where 

m\  and  are  the  mean  and  variance  of  X,  respectively. 
(In  practice,  we  generally  will  not  have  and  at  our 
disposal;  therefore,  their  sample  estimators  must  be  used 
to  estimate  c^.)  Subsequently,  one  may  choose  moment¬ 
matching,  maximum  likelihood  or  so-called  minimum- 
distance  techniques.  It  is  important  to  note  that  maximum 
likelihood  and  minimum-distance  techniques  require  solv¬ 
ing  a  nonlinear  optimization  problem,  whereas  moment¬ 
matching  techniques  require  only  the  true,  or  estimated, 
lower  moments  of  X.  For  this  reason,  we  consider  only 
moment-matching  techniques.  In  what  follows,  let  m/, 
j  =  1, 2,  3,  denote  the  yth  moment  of  X. 

Osogami  and  Harchol-Balter  (2006)  provide  a  nice  sum¬ 
mary  of  the  most  promising  moment -matching  techniques. 


the  majority  of  which  match  the  first  two  or  three  mo¬ 
ments  of  X  to  the  moments  of  a  PH  distribution.  There 
is  little  dispute  (see  Aldous  and  Shepp  (1987))  that  the 
k-phase  Erlang  distribution  provides  very  reliable  approx¬ 
imations  when  <  0.5.  However,  when  >  0.5,  it  is  not 
obvious  which  technique  is  best.  For  example,  Marie  (1980) 
used  two  moments  to  match  a  generalized  Erlang  when 
<  1,  and  matched  two  moments  to  a  C2  distribution 
when  0.5  <  <  1.0.  A  two-moment  matching  algorithm 

using  a  k-phase  Erlang  was  used  for  the  range  0  <  <  0.5. 

Telek  and  Heindl  (2003)  matched  three  moments  of  X  to 
their  two-phase  canonical  acyclic  PH  distribution  when 
0.5  <  <  1.  However,  for  this  range,  their  technique  re¬ 

quires  bounds  on  m3  of  the  form: 

3mj(3c^  —  1  +  V2(l  —  c^)^)  <  m3  <  6m\c^. 

By  contrast,  Marie  (1980)  estimated  the  three  parameters 
of  C2  using  the  first  two  moments,  mi  and  m2,  with 

2  1  1 

/XI  =  — ,  /X2  = - T,  a  =  ^, 

mi  m\c^  2c^ 

resulting  in  a  C2  distribution  representation  with  |3  = 
(1, 0, 0)  and 


Whitt  (1984)  showed  that,  when  >  1,  it  is  important 
to  match  the  first  three  moments  to  reduce  the  maximum 
relative  error  of  the  approximation.  He  used  a  two-branch 
hyper-exponential  distribution,  whereas  Altiok  (1985)  used 
a  C2  distribution,  and  Telek  and  Heindl  (2003)  used  a 
two-phase  canonical  acyclic  PH  distribution  while  match¬ 
ing  three  moments.  The  hyper-exponential  distribution  of 
Whitt  (1984)  requires  the  estimation  of  four  parameters 
while  the  techniques  of  Altiok  (1985)  and  Telek  and  Heindl 
(2003)  require  only  three  parameters.  However,  the  latter 
two  impose  the  requirement: 

3(c2  +  1)2^3 
m3  >  - . 

2 

For  our  semi-Markov  environment  model,  it  is  critical 
to  minimize  the  number  of  parameters  estimated  from  the 
observed  data.  Therefore,  a  PH  distribution  that  uses  a 
minimal  number  of  phases  is  ideal.  To  date,  the  question  of 
determining  a  minimal  representation  with  respect  to  the 
class  of  all  PH  distributions  remains  open.  However,  Os¬ 
ogami  and  Harchol-Balter  (2006)  recently  mapped  general 
distributions  to  representations  that  are  minimal  with  re¬ 
spect  to  the  class  of  acyclic  PH  distributions  (to  which  the 
Ck  and  k-phase  Erlang  distributions  belong)  when  match¬ 
ing  three  moments.  As  previously  noted,  Telek  and  Heindl 
(2003)  examined  canonical  forms,  a  subclass  of  acyclic 
PH  distributions,  that  also  admit  minimal  representations. 
Note  that  the  two-moment  approximation  of  Marie  (1980) 
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Table  1.  PH  approximation  selection  criteria 


Range  of 

PH  approximation 

0  <  c2  <  0.5 

0.5  <  c2  <  1 
c2  >  1 

2-Moment,  k-phase  generalized  Erlang 

2- Moment,  2-phase  Coxian 

3- Moment,  2-phase  Coxian 

for  0  <  <  0.5  provides  a  minimal  representation  when 

matching  only  two  moments. 

By  considering  the  guidance  provided  by  Marie  (1980), 
Whitt  (1984),  Altiok  (1985),  Perros  (1994),  and  Osogami 
and  Harchol-Balter  (2006),  Table  1  summarizes  our  criteria 
for  selecting  the  PH  approximation  type  based  on  the  true, 
or  estimated,  value  of  c^. 

In  what  follows,  we  summarize  the  specific  calculations 
needed  to  match  moments  using  the  three  types  of  PH 
approximations  of  Table  1 .  This  discussion  is  largely  con¬ 
tained  in  Marie  (1980),  Altiok  (1985),  Perros  (1994),  and 
Osogami  and  Harchol-Balter(2006). 


3.2.1.  2-moment,  k-phase  generalized  Erlang  distribution 
When  0  <  <  0.5,  we  will  approximate  the  state  holding 

time  distribution,  H,  by  a  2-moment,  k-phase  generalized 
Erlang  approximation  where  the  integer  k{k>  1)  is  chosen 
such  that  (see  Perros  (1994)): 


1  2  1 

-  <  c  <  - . 

k-  “  (k  -  1) 


Recall  that  a  is  the  probability  that  {0(?)  :  ?  >  0}  tran¬ 
sitions  from  phase  1  to  phase  2,  and  /r  is  the  common 
rate  of  the  k  exponential  phases.  These  two  parameters  are 
obtained  by 


2kc^  +  k-2-{k^  +  A-  Akc'^y^ 
2(c2  +  l)(k  -  1) 

and 


(13) 


1  +  (k  —  \)a 

p  = - 

m\ 


(14) 


where  mi  must  be  estimated  from  observed  data.  Subse¬ 
quently,  a  and  p  are  used  directly  in  the  representation  of 
Equation  (12). 


3.2.2.  2-moment,  2-phase  Coxian  distribution 
When  0.5  <  c2  <  1,  a  2-moment,  2-phase  Coxian  distribu¬ 
tion  is  used  to  approximate  H.  The  representation  Equation 
(10)  requires  the  exponential  rates  pi,  i  =  \,2, . . . ,  k,  and 
the  absorption  probabilities,  (1  -  a,).  Eollowing  a  moment¬ 
matching  algorithm  similar  to  that  for  the  3-moment,  2- 
phase  Coxian  distribution  (discussed  next),  the  parameters 
of  the  2-moment  variant  of  the  2-phase  Coxian  are 

pi=2/mi;  p2  =  \-/m\c~',  a  =  \/2c^, 

where  mi  and  are  estimated  from  the  observed  data.  It  is 
important  to  note  that  using  either  variant  of  the  2-phase 


Coxian  distribution  as  a  PH  approximation  helps  to  limit 
the  growth  of  the  surrogate  environment  state  space. 

3.2.3.  3-moment,  2-phase  Coxian  distribution 

When  c2  >  1,  the  c.d.f  //is  approximated  by  a  3-moment, 

2-phase  Coxian  distribution.  As  noted  by  Whitt  (1984)  and 

Altiok  (1985),  the  third  moment  is  significant  when 

1.  By  taking  successive  derivatives  of  the  EST  of  the  C2 

distribution,  it  has  been  shown  in  Altiok  (1985)  that 

mi^E[Y\  =  —  +  —,  (15) 

Pi  P2 


2i  2(1 -a)  2apip2-2a{pi+p2f 
m2  =  E[Y]  = - 2 - - ’(16) 


p\ 


2  2 

p\pi^ 


and 

m3 

6(1  -  a)  \2apip2{p\  +  P2)  —  (>a{pi  +  p2)^ 


p\ 


3  3 

Mi/^2 


•  (17) 


As  in  Altiok  (1985),  if  we  set  A  =  pi+  /r2and  B  =  pip2, 
substitute,  and  simplify,  the  parameters  of  the  3-moment, 
2-phase  Coxian  distribution  are  given  by 


-  4B 

pi  =  A3 - 2 - ,  p2  =  A-pi, 

a  =  p^^ P2{mipi  -  1).  (18) 


We  statistically  estimate  the  values  mi,  m2,  and  m3  and 
use  these  estimates  in  equations  Equation  (15)  to  (18)  to 
obtain  parameters  pi,  p2,  and  a  for  use  in  representation 
Equation  (10). 

In  the  next  section,  we  present  our  formalized  procedure 
for  approximating  environment  state  holding  times  and  us¬ 
ing  these  within  a  Markovian  framework  to  obtain  lifetime 
distribution  approximations.  Subsequently,  the  quality  of 
these  approximations  will  be  assessed  through  numerical 
experiments  in  Section  5. 


4.  Approximation  procedure 

In  this  section,  we  formalize  our  procedure  for  construct¬ 
ing  a  surrogate  Markovian  environment  process  to  replace 
the  observed  semi-Markov  environment.  To  this  end,  let  us 
clarify  a  few  assumptions.  Eirst,  we  assume  the  availability 
of  sensor  data  that  provides  information  about  the  current 
state  of  the  environment  (e.g.,  load,  speed,  temperature, 
pressure,  humidity,  etc.)  in  which  the  unit  resides,  or  under 
which  it  operates,  while  the  degradation  status  is  unobserv¬ 
able.  Second,  the  environment  is  assumed  to  evolve  as  a 
finite,  irreducible  SMP  that  directly  influences  the  degra¬ 
dation  process  {X{t) :  ?  >  0};  however,  no  assumptions  are 
made  regarding  the  probability  law  of  the  latter.  Third,  we 
assume  that  the  environment’s  true  state  space  can  be  par¬ 
titioned  into  K  distinct  states  such  that  5  =  {1,  2, . . . ,  Ai}, 
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2  <  K  <  oo.  However,  if  K  is  not  known  a  priori,  it  can 
be  estimated  by  using  a  simple  clustering  scheme,  such  as 
that  described  in  Kharoufeh  and  Cox  (2005).  Associated 
with  each  environment  state  is  a  known,  positive  degrada¬ 
tion  rate,  r{j),  so  that  when  Z{t)  =  j,  the  system  degrades 
linearly  at  rate  r{j).  When  Z  can  be  modeled  as  a  CTMC, 
Equations  (6)  and  (7)  require  the  environment’s  initial  dis¬ 
tribution  vector  «,  the  diagonal  matrix  of  degradation  rates 
Rd,  and  the  inhnitesimal  generator  matrix  Q.  However,  be¬ 
cause  we  relax  the  Markovian  assumption  here,  our  aim 
is  to  obtain  a  surrogate  Markovian  environment  Z  with 
state  space  S  and  generator  matrix  Q  by  approximating 
each  state  holding  time  distribution  by  a  2-  or  k-phase  PH 
distribution,  depending  on  the  estimated  value  of  c^,  as  de¬ 
scribed  in  Section  3.  Because  each  PH  approximation  uses 
at  least  two  phases,  we  have  that  1 5|  >  1 5| .  For  instance, 
if  |5|  =  3,  and  each  holding  time  distribution  is  approx¬ 
imated  by  a  2-phase  j’H  distribution,  then  the  surrogate 
environment  process  Q  will  have  3  -h  2  x  3  =  9  states.  Note 
that  the  degradation  rate  matrix,  Rd,  must  also  be  replaced 
byRd,  which  has  the  same  dimension  as  Q. 

Assume  the  homogeneous  environment  process  {Z{t)  : 
?  >  0}  is  perfectly  observable  over  some  time  interval  [0,  r], 
and  that  the  unit  has  not  failed  in  this  interval.  The  en¬ 
vironment  is  continuously  observed  up  to  time  x,  and  at 
each  transition  epoch,  we  record  the  current  and  subse¬ 
quent  state  of  the  random  environment.  Let  R{i,  j)  denote 
the  true  (unknown)  rate  at  which  the  environment  transi¬ 
tions  from  state  i  e  Sto  state  j  e  S,  j  ^  i.  Let  N-di,  j)  be 
the  number  of  transitions  from  i  to  j  in  time  t,  let  A)(/) 
be  the  number  of  visits  to  state  i,  and  let  Wdi)  be  the  total 
time  spent  by  Z in  state  i  during  [0,  t].  It  can  be  shown  (cf 
Basawa  and  Rao  (1980))  that,  for  each  i  e  S  and  j  i,  as 
T  ^  oo: 


NrjiJ) 

Wdi) 


R(iJ)  a.s., 


D(i)  weakly  as  n  ^  oo,  where  D(i)  is  the  unconditional 
holding  time  in  state  i  with  proper  c.d.f  Hi.  For  each  i  e  S, 
the  kth  moment  and  variance  of  D{i)  are  estimated  by 


and 


NAi) 


-j 


n=\ 


NAi) 


1 


(22) 


(23) 


respectively.  Note  that  Equation  (22)  is  the  sample  /rth 
moment  of  the  holding  time  in  state  i  during  [0,  t], 
and  Equation  (23)  is  the  sample  variance  of  the  holding 
time.  The  PH  distribution  selection  procedure  described 
in  Section  3  requires  the  sample  squared  coefficient  of 
variation: 


m\{iy 


i  e  5. 


(24) 


We  now  formalize  the  full  procedure  to  obtain  Z,  Q,  and 

Rd 


Step  0.  Initialization. 

Select  a  sufficiently  large  observation  time  t. 

Step  1.  Observe  the  environment  process  on  [0,  t]. 

During  the  observation  period,  for  each  {i,  j)  such 
that  j  i,  record  the  number  of  transitions  from 
i  j,  Nr(i,  j).  Also  record  the  duration  of  each 
visit  to  state  i,  Dnii)  for  n  =  1,  2, . . . ,  Ndi)-  The 
total  time  spent  in  state  i  during  [0,  t]  is  approxi¬ 
mated  by 


Rd) 

Wrii)  ^ 

n=l 


Therefore,  for  r  sufficiently  large,  we  approximate  R(i,  j) 
by 

R(i,  j)  Rr{i,  j)  =  j  +  i,  (19) 

Writ) 

and  set 

Rr(i,  0  =  -  I]  Mi,  j),  i  e  (20) 

where  Rr(i,  i)  <0.  Now  let  P  denote  the  estimated  transi¬ 
tion  matrix  of  {£„  :  n  >  0},  the  Markov  chain  embedded  at 
environment  transitioi^epochs.  The  statistical  estimator  of 
the  (i,  7)th  element  of  P  is 


'i.j 


RAi,j) 

-RAi.i) 

0 


if  j  7^  i, 
if  j  =  i. 


(21) 


For  n  >  1  and  i  e  S,  let  D„(/)  be  the  duration  of  the 
environment’s  nth  visit  to  state  i,  and  assume  that  Dnii) 


since  W-di)  is  (almost  surely)  bounded  above  by  the 
sum  on  the  right-hand  side. 

Step  2.  Estimate  the  required  parameters. 

Using  Ndi),  Dnii),  n  =  1, . . . ,  Nd^Ndi,  j)  and 
Wdi),  estimate  the  transition  matrix  P  using  Equa¬ 
tions  (19)  to  (21).  For  each  i  e  S,  compute  the  first 
three  sample  moments  nikii),  k=  1,  2,  3,  the  sam¬ 
ple  variance  sf,  and  the  sample  squared  coefficient 
of  variation  using  Equations  (22)  to  (24). 

Step  3.  Select  the  PH  distribution  approximation  for  Hi. 

Using  nikii),  k=  1 ,  2,  3,  and  c],  i  e  S,  select  the  ap¬ 
propriate  PH  distribution  approximation  in  accor¬ 
dance  with  Table  1.  Using  the  moment-matching 
methods  described  in  Section  3,  compute  the  ma¬ 
trix  T,  for  each  i  e  S.  Note  that  T®  follows  directly 
from  T,  since  T,e  +  T®  =  0.  Therefore,  for  each 
t  >  0,  the  PH  approximation  of  Hi  is  given  by 

Hiit)  =  1  -  |3  exp(T,  Oe,  i  e  S. 
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Step  4.  Construct  surrogate  environment  process. 

Suppose  1 5|  =  Ai  and  that  Hi  is  approximated  by  a 
PH  distribution  with  ki  phases.  That  is,  an  absorb¬ 
ing  Markov  chain  with  ki  +  1  states  replaces  the 
original  state  i.  Th^efore,  the  surrogate  Marko¬ 
vian  environment,  {Z{t)  :  ?  >  0},  has  state  space  S 
with  cardinality: 

K 

\S\  =  K  +  Y,ki. 

i=\ 

For  this  reason,  it  is  important  to  use  the  smallest 
possible  integer  k,  for  each  i  e  5^Next,  construct 
the  surrogate,  5- valued  CTMC  {Z{t)  :  ?  >  0}  with 
generator  matrix  Q  given  by 


f  An 

Ai2 

•  Ai;,\ 

A21 

A22 

A2^: 

Q  = 

vA^:i 

A^:2 

•  Skk  ) 

where 

/T, 

T°  \ 

A,','  = 

c  Rr(i,  i)  j 

1“ 

for  j 

i  •"  I 

\cRdi,j)  oy 


and  c  is  a  very  large,  positive  real  number.  The 
parameter  c  is  needed  in  the  numerical  implemen¬ 
tation  to  ensure  that  {Z{t)  :  t  >  0}  instantaneously 
transitions  to  the  next  state  chosen  by  P  when  the 
environment’s  sojourn  in  state  /  e  5  is  complete. 
Note  that  T,  is  a  k,  x  k,  matrix  while  T°  is  a  k,  - 
dimensional  column  vector.  The  degradation  rate 
matrix  must  also  be  expanded  accordingly.  The  new 


matrix  is  of  the  form 

/A(r(l)) 

0 

0  \ 

0 

A(r(2))  .. 

0 

Rd  = 

1  0 

0 

where  A(a)  is  a  diagonal  matrix  whose  diagonal 
entries  are  all  a. 

Step  5.  Approximate  the  lifetime  distribution  function  and 
moments. 

By  replacing  the  SMP  environment  {Z{t)  :  ?  >  0} 
with  state  spa(x  S  and  generator  Q  by  the  CTMC 
environment  {Z{t)  :  ?  >  0}  with  state  space  S  and 
generator  Q,  the  LST  of  the  unconditional  lifetime 
distribution  is  given  by 

F{u,  0=1-  aexp[(Q  -  wRdjtje,  (25) 


and  its  conditional  counterpart  is 

Fi{u,  t)  =  l  -  e'  exp[(Q  -  MRd)?]e.  (26) 

Moreover,  the  LST  of  the  nth  moment  of  77,  con¬ 
ditioned  on  the  initial  state  of  the  environment,  is 
obtained  by  (see  Kharoufeh  et  al.  (2006)) 

H[t:]  =  file'  (uRd  -  Qf"  e,  n  >  1.  (27) 


While  the  procedure  described  here  appears  somewhat 
involved,  it  has  the  advantage  of  requiring  only  event  times 
and  count  data  for  the  environment  process.  There  is  no 
need  for  failure  time  observations,  or  for  degradation  ob¬ 
servations  (besides  those  that  will  be  needed  initially  to 
estimate  the  rates  r(J),  j  e  S).  The  next  section  provides 
two  numerical  experiments  to  assess  the  quality  of  our  ap¬ 
proximations  and  to  illustrate  the  viability  of  the  approach. 


5.  Numerical  experiments 

This  section  provides  two  extensive  numerical  experiments 
to  assess  the  quality  of  distribution  approximations  ob¬ 
tained  by  the  procedures  described  in  Sections  3  and  4. 
In  lieu  of  real  data  describing  the  evolution  of  a  random 
environment,  we  randomly  generated  1000  test  scenarios 
for  each  of  the  two  experiments.  Specifically,  we:  (i)  chose 
holding  time  distributions  (H)  for  SMP  environments  and 
randomized  their  parameter  values;  (ii)  randomly  gener¬ 
ated  the  transition  probability  matrix  (P)  of  the  discrete¬ 
time  Markov  chain  embedded  at  transition  epochs;  and  (iii) 
randomly  generated  the  degradation  rate  (r(i))  associated 
with  the  fth  environment  state.  To  the  outside  observer,  the 
SMP  process  that  drives  the  environment’s  evolution  is  un¬ 
known.  Therefore,  we  observed  only  the  states  visited  by 
the  environment  and  the  time  spent  in  each  state  over  a 
time  interval  [0,  t].  Subsequently,  the  statistical  estimators 
described  in  Section  4  were  calculated.  From  these  values, 
PH  distribution  approximations  and  the  surrogate  (Marko¬ 
vian)  environment  process  were  constructed  to  approximate 
the  c.dj".  of  77. 

Let  Fi(x,  t)  denote  the  approximate  lifetime  c.d.f  gener¬ 
ated  by  the  surrogate  environment  [Z{t)  :  t^  0}  with  gen¬ 
erator  matrix  Q  and  degradation  matrix  Rd.  Denote  by 
Fi{x,  t)  the  c.d.f  of  77  obtained  by  simulating  the  degrada¬ 
tion  process,  {X{t) :  ?  >  0},  until  it  first  reaches  the  critical 
threshold  value  x.  For  each  of  the  1000  scenarios,  and  for 
both  experiments,  77  (x,  t)  was  constructed  using  20  000  ob¬ 
servations  of  77,  given  that  Z(0)  =  i.  Therefore,  we  view  the 
simulated  c.d.f  as  the  ccmpletely  specified,  hypothesized 
benchmark.  To  compare  77  (x,  t)  and  77  (x,  t),  we  employed 
the  non-parametric  two-sided  Kolmogorov-Smirnov  (KS) 
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goodness-of-fit  test,  which  tests  the  hypothesis: 

Hq  :  Fi{x,  t)  =  Fi(x,  t),  for  all  t  e  R, 

Hi  :  Fi{x,  t)  7^  Fi{x,  t),  for  at  least  one  t  e  R. 

The  distribution  T’,  (x,  t)  was  sampled  on  a  finite  set  T  c 
=  [0,  oo)  whose  range  depends  on  observed  simulated 
outcomes.  The  test  statistic  for  the  two-sided  KS  test  (see 
Conover  (1999),  p.  431).  is  given  by 

8  =  sup  \  Fi{x,  t)  -  Fi(x,  ?)|  , 

the  maximum  absolute  deviation  between  the  two  distribu¬ 
tions  over  the  set  T.  The  null  hypothesis,  Hq,  is  rejected  at 
the  a  level  of  significance  if  £  >  e*  where  £*  is  the  (1  —  a) 
quantile  of  the  test  statistic.  Let  n  =  |  r|  denote  the  cardi¬ 
nality  of  T.  At  the  0.05  level,  the  critical  value  £*  is  given 
by 

£*  =  1.36n^^/^ 

if  n  >  40,  and  if  n  <  40,  the  critical  value  is  obtained  from 
Conover  (1999,  Table  A1 3,  p.  547).  All  hypothesis  tests  were 
performed  at  the  a  =  0^5  level  of  significance. 

Numerical  values  of  (x,  t)  were  obtained  via  the  inverse 
Laplace  transform: 

Fiix,t)  =  R~\u-^Fi{u,t)).  (28) 

where  ij  (M,  t)  is  given  by  Equation  (26),  and  denotes 
the  inverse  Laplace  transform  operator.  For  illustrative  pur¬ 
poses,  we  also  assess  the  quality  of  the  lower  moment  ap¬ 
proximations  given  by 

E[r;|Z(0)  =  /]  =  [u-^Ut:])  ,  «  >  l,  (29) 

where  E,  [TjJ*]  is  given  by  Equation  (27).  The  inverse  Eaplace 
transforms  (28)  and  (29)  were  obtained  by  coding  a  numer¬ 
ical  Eaplace  transform  inversion  algorithm  due  to  Abate 
and  Whitt  (1995).  This  algorithm  and  the  discrete-event 
simulation  models  were  coded  in  the  MATEAB®  com¬ 
puting  environment  and  executed  on  a  personal  computer 
equipped  with  an  Intel®  Core^“  2  Duo  CPU  operating  at 
3.00  GHz  with  2.00  GB  of  RAM. 

5.1.  Experiment  1:  degradation  of  a  turbine  blade 

Consider  an  aircraft  engine  turbine  blade  that  operates  at 
high  temperatures  and  experiences  centrifugal  stresses.  The 
rotational  speed  of  the  blade  varies  depending  on  uncer¬ 
tain  prevailing  conditions.  For  example,  the  engine  load 
varies  as  the  aircraft  is  exposed  to  different  flight  condi¬ 
tions  (i.e.,  takeoff,  maximum  climb,  maximum  cruise,  loi¬ 
ter,  flight  idle,  taxi,  ground  idle  and  cutoff),  or  at  different 
flight  altitudes.  Moreover,  if  inclement  weather  is  encoun¬ 
tered  mid-flight,  the  aircraft  may  change  altitude,  acceler¬ 
ate,  or  decelerate  to  minimize  the  effects  of  turbulence.  The 
induced  centrifugal  stresses  can  lead  to  elongated  particles 
at  the  microstructure  level  that  can  result  in  degradation  of 
fatigue  strength  that  results  in  voids  and  crack  initiation  in 
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Table  2.  Summary  input  data  for  the  turbine  blade  experiment 


State 

Rotation 
speed  ( rpm ) 

Holding  time 
distribution 

Degradation 
growth  rate 

1 

1000-4999 

Beta  (a,  b) 

r(l)  ~  t/(0,  5) 

2 

5000-8999 

Beta  (a,  b) 

r(2)  ~  t/(0,  5) 

3 

9000-12999 

Weibull  (c,  d) 

r(3)  ~  G(0,  5) 

4 

13000-15999 

Weibull  (c,  d) 

r(4)  ~  t/(0,  5) 

the  blade  and  continued  loading  exacerbates  the  degrada¬ 
tion.  This  phenomenon  has  been  documented,  for  exam¬ 
ple,  in  Persson  and  Persson  (1993),  Tamarin  (2002),  Ejaz 
and  Tauqir  (2006).  The  blade’s  dynamic  and  stochastic  ro¬ 
tational  speed  can  be  viewed  as  a  proxy  for  the  induced 
stresses.  Suppose  the  degradation  process  {U(?) :  ?  >  0} 
tracks  the  degradation  level  (e.g.,  the  length  of  a  crack) 
of  the  blade,  and  let  Z{t)  denote  the  rotational  speed  of  the 
blade  at  time  t.  We  partition  the  possible  range  of  speeds 
into  four  non-overlapping  intervals  that  describe  four  dis¬ 
tinct  operating  states  (or  flight  conditions).  The  state  space 
of  the  environment  is  5=  {1,  2,  3,  4}.  The  critical  thresh¬ 
old  for  the  degradation  level  is  x  =  20.0  units.  The  state 
descriptions,  holding  time  distributions,  and  degradation 
growth  rates  are  summarized  in  Table  2. 

In  Table  2,  Beta  (a,  b)  denotes  a  beta  probability  distri¬ 
bution  with  shape  parameters  a  and  Z).  (Note  that  the  beta 
distribution  does  not  have  a  closed-form  c.d.f).  The  pa¬ 
rameter  values  a  and  b  were  randomly  generated  for  each 
of  the  1000  scenarios  so  that: 

a~U(L5)  and  b-^U{\,5). 

This  range  of  values  allows  for  a  variety  of  distribution 
shapes  including  right-skewed,  left-skewed  and  symmetric 
forms.  Weibull  (c,  d)  denotes  a  Weibull  distribution  with 
shape  parameter  c  and  scale  parameter  d.  These  parameter 
values  were  also  randomly  generated  for  each  of  the  1000 
scenarios  as  follows: 

c  ~  U(0.5,  6)  and  d  ~  U(0.5,  6). 

The  transition  probability  matrix  of  the  SMP  environment 
was  generated  by  first  creating  a  4  x  4  matrix  B  =  [bij]  such 
that  bij  ~  U(0,  1),  j  7^  i,  and  =  0,  which  was  subse¬ 
quently  rescaled  to  ensure  a  stochastic  matrix.  Specifically, 
the  (/,  7)th  element  of  the  randomly  generated  transition 
matrix  P  is  obtained  by 

PiJ  =  ^4  ,  ’  07  e  5, 

for  each  of  the  1000  scenarios.  From  Table  2,  it  is  seen  that 
the  matrix  of  degradation  rates  is  also  randomly  generated 
in  each  case. 

We  assume  Z  begins  in  the  first  state  with  a  proba¬ 
bility  of  one;  therefore,  the  initial  probability  vector  is 
ej  =  (1, 0,  0, 0).  We  simulated  Z and  observed  its  evolution 
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Fig.  1.  Simulated  versus  approximated  lifetime  c.d.f:  turbine  blade  experiment. 


for  T  =  10  000  time  units.  From  the  simulated  environment, 
we  estimated  the  parameters  P,  mkii),  sf,  cj,  i  =  1,2,3,  4, 
k=  1, 2,  3,  using  Equations  (19)  to  (24).  Subsequently,  we 
applied  the  criteria  of  Table  Ho  select  PH  distributions  for 
each  state  i  and  constructed  Q  and  Ra  using  the  procedure 
described  in  Section  4.  For  the  sake  of  brevity,  we  do  not 
include  here  the  PH  representations  (|3,  T,),  /  =  1,  . . . ,  4. 

In  988  of  1000  randomized  scenarios,  we  fail  to  reject  the 
null  hypothesis  that  Fi{x,  t)  is  equivalent  to  the  benchmark 
c.d.f  Fi{x,  t)  at  the  a  =  0.05  level.  Of  the  12  cases  in  which 
we  reject  i/o,  six  were  the  result  of  instability  in  the  numer¬ 
ical  Laplace  transform  inversion  algorithm.  Specifically,  in 
the  rare  cases  when  the  c.d.f  is  not  sufihciently  smooth,  the 
inversion  algorithm  of  Abate  and  Whitt  (1995)  may  not  per¬ 
form  well  without  modifications.  Of  the  remaining  six  cases 
for  which  //q  was  rejected,  the  average  maximum  absolute 
deviation  in  probability  was  approximately  0. 1478.  In  sum¬ 
mary,  we  fail  to  reject  nearly  99%  of  the  1000  randomized 
scenarios. 

As  a  sample  illustration  of  the  quality  of  the  approxi¬ 
mations,  Fig.  1  depicts  the  benchmark  (simulated)  lifetime 
c.d.f,  Fi{x,  t),  and  the  approximated  lifetime  c.d.f  Fj^x,  t) 
obtained  using  the  surrogate  Markov  environment  Q  and 
degradation  rate  matrix  Rd.  In  this  particular  scenario, 
our  procedure  yielded  a  41  x  41  generator  matrix  (Q).  The 
maximum  absolute  deviation  in  probability  of  this  scenario 


is  e  0.02562.  These  results  indicate  that  the  surrogate 
environment  process  with  PH-approximated  holding  times 
very  closely  tracks  the  simulated  lifetime  c.d.f,  even  when 
holding  times  are  clearly  non-memoryless  and  their  param¬ 
eters  are  randomized. 

We  also  compared  the  first  and  second  moments  of  the 
full  lifetime  distribution,  given  that  the  environment  starts 
in  state  1.  These  values  are  summarized  in  Table  3. 

The  results  of  Table  3  illustrate  the  quality  of  the  lower 
moment  approximations  obtained  using  our  procedure. 
Both  values  differ  from  their  simulated  counterparts  by 
less  than  1%. 

5.2.  Experiment  2:  chemical  coating  decomposition 

Now  consider  a  protective,  automotive  coating  that  is  sub¬ 
ject  to  weathering  by  the  outdoor  environment.  Environ¬ 
mental  effects,  such  as  temperature  and  solar  radiation. 

Table  3.  Lifetime  moment  comparisons  for  experiment  1 

PFl  Percentage 

Parameter  Simulated  approximation  error 

E[F\Z{G)=\]  4.763766782  4.777832971  0.29527453 

E[T^\Z{Q)=\]  22.77583335  22.91001248  0.589129367 
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Table  4.  Summary  input  data  for  the  chemical  coating  experiment 


State 

Sky 

condition 

Temperature 

Holding  time 
distribution 

Decomposition 

rate 

1 

Cloudy 

<32°F 

Weibull  (c,  d) 

r(l) 

~  [7(0,  2) 

2 

Cloudy 

>32°F 

Beta  (a,  b) 

r(2) 

~  67(0,  2) 

3 

Sunny 

<32°F 

Weibull  (c,  d) 

r(3) 

~  67(0,  2) 

4 

Sunny 

>32°F 

Beta  (a,  b) 

r(4) 

~  67(0,  2) 

5 

Rain 

>32°F 

Gamma  {g,h) 

r(5) 

~  67(0,  2) 

cause  chemical  decomposition  (degradation)  of  the  coat¬ 
ing  that  can  be  measured  in  terms  of  gloss  loss  and/or 
color  change.  The  failure  time  can  be  defined  as  the  first 
time  that  the  cumulative  degradation  of  the  coating  reaches 
a  critical  threshold.  Here,  we  assume  the  rate  of  degrada¬ 
tion  of  the  coating  depends  on  the  ambient  environment 
to  which  it  is  exposed,  and  the  critical  threshold  is  x  =  5.0 
units.  The  environment  is  characterized  by  hve-state  semi- 
Markov  processes  that  are  defined  in  Table  4  along  with  the 
degradation  rates  and  state  holding  time  distributions. 

In  Table  4,  Beta  {a,  h)  denotes  a  beta  probability  distri¬ 
bution  with  shape  parameters  a  and  h  and  Weibull  (c,  d) 
denotes  a  Weibull  distribution  with  shape  parameter  c  and 
scale  parameter  d.  Gamma  {g,h)  denotes  a  gamma  dis¬ 
tribution  with  parameters  g  and  h.  For  each  of  the  1000 
scenarios  the  distribution  parameter  values  were  randomly 


generated  as  follows: 

a,h--  t/(0.1,  3.1), 
c,i  ~  U{1,  3), 
g,h  ~  67(0.1,2.1), 

and  the  5x5  transition  probability  matrices  of  the  SMP 
environments  were  generated  randomly  as  in  the  turbine 
blade  experiment.  Here,  we  again  assume  that  Z  begins  in 
state  1  (with  a  probability  of  one)  so  that  the  initial  probabil¬ 
ity  vector  is  e)  =(1,0,  0, 0).  We  simulated  Z  and  observed 
its  evolution  for  t  =  10  000  time  units.  From  the  simulated 
environment,  we  estimated  the  parameters  P,  sf,  c^, 
i  =  1,  2,  3,  4,  5,  A:  =  1,  2,  3,  using  Equations  (19)  to  (24). 
We  used  cf  to  determine  the  appropriate  PH  distribution 
^proximation  for  each  7  e  Susing  Table  1  and  constructed 
Q  and  Rd  using  the  algorithm  described  in  Section  4. 

Of  the  1000  randc^ized  scenarios,  we  fail  to  reject  the 
null  hypothesis  that  Fi(x,  t)  is  equivalent  to  F\{x,  t)  in  980 
cases  at  the  a  =  0.05  level.  Of  the  20  cases  for  which  we 
reject  //q,  ten  were  caused  by  instability  in  the  numeri¬ 
cal  Laplace  transform  inversion  algorithm  as  in  the  hrst 
experiment.  For  the  remaining  ten  cases,  the  average  max¬ 
imum  absolute  deviation  in  probability  was  approximately 
0.2230.  To  summarize  the  results  of  this  experiment,  we  fail 
to  reject  the  null  hypothesis  that  our  approximated  c.d.f  is 


Fig.  2.  Simulated  versus  approximated  lifetime  c.d.f:  chemical  coating  experiment. 


Semi- Markov  models 
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Table  5.  Lifetime  moment  comparisons  for  experiment  2 


PH 

Percentage 

Parameter 

Simulated 

approximation 

error 

E[TAZ{(S)  = 

1] 

4.290965585 

4.249862484 

0.957898641 

E[Tl\Z(G)  = 

1] 

19.0983667 

18.74008389 

1.875986669 

equivalent  to  the  benchmark  c.d.f.  in  98%  of  the  random¬ 
ized  scenarios. 

To  further  illustrate  the  quality  of  the  approximations, 
Fig.  2  depicts  the  simulated  hfetime  c.d.f,  Fi{x,  t),  and  the 
approximated  lifetime  c.d.f  F^{x,  t)  obtained  using  the  sur¬ 
rogate  Markov  environment  Q  and  degradation  rate  matrix 
Rd.  For  this  particular  experiment,  our  procedure  yielded 
a  generator  matrix  Q,  which  is  82  x  82,  or  twice  as  large  as 
the  illustrative  case  in  experiment  1.  Nevertheless,  the  re¬ 
sulting  maximum  absolute  deviation  in  probability  of  this 
example  is  e  0.02648. 

The  lower  two  moments  of  the  hfetime  distribution  were 
also  compared,  and  the  results  are  summarized  in  Table  5. 

As  for  the  first  example,  the  approximation  of  the  lower 
moments  is  also  of  very  high  quality.  The  maximum  percent 
difference  is  under  2%. 

While  the  results  of  these  two  randomized  experiments 
are  promising,  it  will  be  instructive  to  assess  the  quality  of 
our  approximations  using  real  data.  Unfortunately,  such 
data  were  unavailable  for  use  in  this  study.  We  further  elab¬ 
orate  upon  the  advantages  and  limitations  of  our  approach 
in  Section  6. 


6.  Conclusions 

In  this  article,  we  have  presented  a  novel  technique 
for  incorporating  environmental  effects  on  component 
degradation  for  the  purpose  of  evaluating  reliability  in¬ 
dices.  The  modeling  framework  and  associated  numeri¬ 
cal  techniques  provide  improved  flexibility  by  allowing  for 
semi-Markovian  environments  and  general  degradation 
patterns.  These  characteristics  make  the  approach  appeal¬ 
ing  whenever  it  is  possible  to  assess  physical  degradation  as 
a  function  of  the  temporal  evolution  of  the  environment. 
Within  our  framework,  the  environment  can  assume  a  va¬ 
riety  of  roles  (e.g.,  time-varying  operating  conditions,  the 
ambient  environment,  etc.). 

While  the  procedures  described  in  Sections  3  and  4  are 
mathematically  sound  and  relatively  easy  to  implement, 
they  do  impose  moderately  restrictive  assumptions,  namely 
that:  (i)  the  environment  evolves  in  a  temporally  homo¬ 
geneous  manner;  (ii)  the  number  of  distinct  states  {K)  is 
known;  (iii)  the  environment  transitions  between  states  ac¬ 
cording  to  a  Markov  chain;  and  (iv)  the  future  degradation 
of  the  unit  is  independent  of  the  history  of  the  degrada¬ 
tion.  Despite  these  limitations,  the  approximations  provide 
an  important  extension  to  prior  methods  that  either  ignore 


the  effects  of  the  environment  entirely,  or  assume  that  the 
environment  is  completely  memoryless. 

In  the  future,  it  will  be  instructive  to  consider  models 
for  systems  whose  rates  of  degradation  depend  not  only  on 
the  state  of  the  environment  but  also  on  the  current  level 
of  degradation.  It  may  also  be  useful  to  consider  Bayesian 
techniques  for  updating  the  parameters  of  the  environment 
as  additional  sensor  data  becomes  available.  If  the  pro¬ 
posed  models  are  to  be  of  any  practical  value  to  engineers, 
it  is  vital  to  provide  an  accurate  estimate  of  the  mapping  r 
that  describes  the  evolution  of  degradation  as  a  function  of 
the  environment.  To  this  end,  real  degradation  data  are  re¬ 
quired,  as  is  the  guidance  and  experience  of  subject  matter 
experts,  to  ensure  the  degradation  rates  are  selected  in  ac¬ 
cordance  with  known  physical  principles.  Finally,  as  noted 
in  Section  5,  it  will  be  important  to  evaluate  the  procedure 
using  real  environmental  and  degradation  data. 
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